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I will review the progress toward a finite baryon density algorithm in the 
canonical ensemble approach which entails particle number projection from 
the fermion determinant. These include an efficient Pade-Z2 stochastic esti- 
mator of the Tr log of the fermion matrix and a Noisy Monte Carlo update 
to accommodate unbiased estimate of the probability. Finally, I will propose 
a Hybrid Noisy Monte Carlo algorithm to reduce the large fluctuation in the 
estimated Tr log due to the gauge field which should improve the acceptance 
rate. Other application such as treating u and d as two separate flavors is 
discussed. 

1 Introduction 

Finite density is a subject of keen interest in a variety of topics, such as 
nuclear equation of state, neutron stars, quark gluon plasma and color su- 
perconductivity in nuclear physics and astrophysics, and high temperature 
superconductors in condensed matter physics. Despite recent advance with 
small chemical potential at finite temperature yQ, the grand canonical ap- 
proach with chemical potential remains a problem for finite density at zero 
temperature. 

The difficulty with the finite chemical potential in lattice QCD stems from 
the infamous sign problem which impedes important sampling with positive 
probability. The partition function for the grand canonical ensemble is repre- 
sented by the Euclidean path-integral 



where the fermion fields with fermion matrix M has been integrated to give 
the determinant. U is the gauge link variable and S g is the gauge action. The 
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chemical potential is introduced to the quark action with the e^ a factor in 
the time-forward hopping term and e~^ a in the time-backward hopping term. 
Here a is the lattice spacing. However, this causes the fermion action to be 
non-75-Hermitian, i.e. 75M75 ^ M>. As a result, the fermion determinant 
detM[U] is complex that leads to the sign problem. 

2 Finite Chemical Potential 

There are several approaches to avoid the sign problem. It was proposed by 
the Glasgow group 2 that the sign problem can be circumvented based on the 
expansion of the grand canonical partition function in powers of the fugacity 
variable e M ' T , 

B=3V 

Z GC (fi/T,T,V) = e^ TB Z B (T,V), (2) 

B = ~3V 

where Zb is the canonical partition function for the baryon sector with baryon 
number B. Zqq is calculated with reweighting of the fermion determinant 
Since Zqc(/J>/T,T, V) is calculated with reweighting based on the gauge con- 
figuration with /i = 0, it avoids the sign problem. However, this does not work, 
except perhaps for small /i near the finite temperature phase transition. We 
will dwell on this later in Sec. |3| This is caused by the 'overlap problem' |3] 
where the important samples of configurations in the fi = simulation has ex- 
ponentially small overlap with those relevant for the finite density. To alleviate 
the overlap problem, a reweighting in multi-parameter space is proposed 0] 
and has been applied to study the end point in the T-fi phase diagram. In 
this case, the Monte Carlo simulation is carried out where the parameters in 
the set do include fi = and /3 C which corresponds to the phase transition 
at temperature T c . The parameter set a in the reweighted measure include 
mu 7^ and an adjusted p in the gauge action. The new /? is determined from 
the Lee- Yang zeros so that one is following the transition line in the T-/x plane 
and the large change in the determinant ratio in the reweighting is compen- 
sated by the change in the gauge action to ensure reasonable overlap. This 
is shown to work to locate the transition line from /i = and T — T c down 
to the critical point on the 4 4 and 6 3 x 4 lattices with staggered fcrmions 0]. 
While the multi-parameter reweighting is successful near the transition line, 
it is not clear how to extend it beyond this region, particularly the T = case 
where one wants to keep the j3 and quark mass fixed while changing the \i. 
One still expects to face the overlap problem in the latter case. It is shown jS] 
that Taylor expanding the observables and the rewriting factor leads to coef- 
ficients expressed in local operators and thus admits study of larger volumes, 
albeit still with small fi at finite temperature. 

In the imaginary chemical potential approach, the fermion determinant 
is real and one can avoid the sign problem |SJ EJ. In practice, a ref- 
erence imaginary chemical potentials is used to carry out the Monte Carlo 
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calculation and the determinants at other chemical potential values are cal- 
culated through a bosonic Monte Carlo calculation so that one can obtain 
the finite baryon partition function Zb(T,V) through the Fourier transform 
of the grand canonical partition function Zccin/T, T,V) [HI- However, this is 
problematic for large systems when the determinant cannot be directly calcu- 
lated and it still suffers from the overlap problem. The QCD phase diagram 
has been studied with physical observables Taylor expanded and analytically 
continued to the real \i jH]. Again, due to the overlap problem, one is limited 
to small real fi near the finite temperature phase transition. 

3 Finite Baryon Density - A Canonical Ensemble 
Approach 

An algorithm based on the canonical ensemble approach to overcome the 
overlap problem at zero temperature is proposed ^0!- To avoid the overlap 
problem, one needs to lock in a definite nonzero baryon sector so that the 
exponentially large contamination from the zero-baryon sector is excluded. 
To see this, we first note that the fermion determinant is a superposition of 
multiple quark loops of all sizes and shapes. This can be easily seen from the 
property of the determinant 



Upon a hopping expansion of logAf, Tr log M represents a sum of single 
loops with all sizes and shapes. The determinant is then the sum of all 
multiple loops. The fermion loops can be separated into two classes. One 
is those which do not go across the time boundary and represent virtual 
quark-antiquark pairs; the other includes those which wraps around the time 
boundary which represent external quarks and antiquarks. The configuration 
with a baryon number one which contains three quark loops wrapping around 
the time boundary will have an energy Mb higher than that with zero baryon 
number. Thus, it is weighted with the probability e - M sN t a t com p arec i w ith 
the one with no net baryons. We see from the above discussion that the fermion 
determinant contains a superposition of sectors of all baryon numbers, pos- 
itive, negative and zero. At zero temperature where MsNtat 3> 1, the zero 
baryon sector dominates and all the other baryon sectors are exponentially 
suppressed. It is obvious that to avoid the overlap problem, one needs to se- 
lect a definite nonzero baryon number sector and stay in it throughout the 
Markov chain of updating gauge configurations. To select a particular baryon 
sector from the determinant can be achieved by the following procedure ^5 : 
first, assign an U(i) phase factor e~ % ^ to the links between the time slices t 
and t + 1 so that the link U/U^ is multiplied by e~ 1 ^ je % ^\ then the particle 
number projection can be carried out through the Fourier transformation of 
the fermion determinant like in the BCS theory 




(3) 



n=l 
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Pn = ^~ d<pe- i<t,N 'detM[(j)} (4) 
2tt Jo 

where N is the net quark number, i.e. quark number minus antiquark number. 
Note that all the virtual quark loops which do not reach the time boundary 
will have a net phase factor of unity; only those with a net N quark loops across 
the time boundary will have a phase factor e*^ which can contribute to the 
integral in Eq. Since QCD in the canonical formulation does not break 
Z(3) symmetry, it is essential to take care that the ensemble is canonical with 
respect to triality. To this end, we shall consider the triality projection |11II12| 
to the zero triality sector 

det M=^ detM[<j> + k2n/3}. (5) 

fc=0,±l 

This amounts to limiting the quark number N to a multiple of 3. Thus the 
triality zero sector corresponds to baryon sectors with integral baryon num- 
bers. 

8 3 x 12 (3 = 6.0 k = 0.150 Wilson action lattice 



j i i i i 

2 4 6 8 

<j> = 2ir * x/9 

Fig. 1. Tr log M [4>] for a 8 3 x 12 configuration with Wilson action as a function of 

<j>. 

Another essential ingredient to circumvent the overlap problem is to stay 
in the chosen nonzero baryon sector so as to avoid mixing with the zero baryon 
sector with exponentially large weight. This can be achieved by preforming 
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(j> = 2ir * x/9 

for a 8 x 20 2 x 4 finite temperature configuration with dynamical 



the baryon number projection as described above before the accept/reject step 
in the Monte Carlo updating of the gauge configuration. If this is not done, 
the accepted gauge configuration will be biased toward the zero baryon sector 
and it is very difficult to project out the nonzero baryon sector afterwords. 
This is analogous to the situation in the nuclear many-body theory where it 
is known J3| that the variation after projection (Zeh-Rouhaninejad-Yoccoz 
method |141 115) ) is superior than the variation before projection (Peierls- 
Yoccoz method QUI)- The former gives the correct nuclear mass in the case 
of translation and yields much improved wave functions in mildly deformed 
nuclei than the latter. 

To illustrate the overlap problem, we plot in Fig. 1 Tr log M [</>] for a config- 
uration of the 8 3 x 12 lattice with the Wilson action at (3 — 6.0 and n = 0.150 
which is obtained with 500 Z2 noises. We see that the it is rather flat in 4> 
indicating that the Fourier transform in Eq. (0J) will mainly favor the zero 
baryon sector. On the other hand, at finite temperature, it is relatively easier 
for the quarks to be excited so that the zero baryon sector does not neces- 
sarily dominate other baryon sectors. Another way of seeing this is that the 
relative weighting factor e - M BN t a t can ^ e o(l) at finite temperature. Thus, 
it should be easier to project out the nonzero baryon sector from the determi- 
nant. We plot in Fig. 2 a similarly obtained Tr log M[4>] for configuration of 
the 8 x 20 2 x 4 lattice at finite temperature with (3 — 4.9 and k = 0.182. We see 
from the figure that there is quite a bit of wiggling in this case as compared 
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to that in Fig. 1. This implies that it is easier to project out a nonzero baryon 
sector through the Fourier transform at finite temperature. 



4 Noisy Monte Carlo with Fermion Determinant 

In order to implement the canonical ensemble approach, it is clear that one 
needs to evaluate the fermion determinant for the purpose of particle projec- 
tion. Since the pseudofermion approach does not give the determinant in the 
Markov process, it is not applicable. In view of the fact that it is impractical 
to calculate the determinant directly for realistic volumes, a Monte Carlo al- 
gorithm which accommodates an unbiased estimate of the probability and an 
efficient way to estimate the determinant are necessary for the finite baryon 
density calculation. 

A noisy Monte Carlo algorithm [2j with Pade-Z2 estimates ^HIE] of the 
Tr log of the fermion matrix are developed toward this goal and a numeri- 
cal simulation with Wilson dynamical fermion is carried out [201 • We shall 
summarize the progress made so far. 

The QCD partition function can be written in the form 



dUe~ s ^ / P^rn) 
J i=i 

oo 

x n fa pp (^) f( u >v,p) , (6) 

k=2 

where S g {U) is the gauge action. f(U,rj,p) stands for f(U, {fy}, {pk}) which 
is an unbiased stochastic estimator [23 of the fermion determinant e Tr log M 
via an infinite number of auxiliary variables p k and r)i. P v (rii) — S(\rji\ — 1) 
is the distribution for the Z 2 noise rji and P p {pk) = 0(pk) — @{Pk — 1) is the 
flat distribution for < pk < 1. With f(U, {pk}) being the stochastic 
expansion 

f(U, {77,}, { P k}) = 1 + I xi + 6 (1 - p2) jxa + Q - p 3 ) {x 3 + ■ ■ ■ 

\-Pr)&n + ■■•}}}} (7) 

where Xi = rj\ In M(U)rji, one can verify |21| that 

oo oo 

Y[d Vi P^ Vi ) Ud Pk PP(p k )(f(U,{ m },{p k })) = e Tr * M M, ( 8 ) 

i=l k=2 

and the stochastic series terminates after e terms on the average. 



A Finite Baryon Density Algorithm 



7 



Since the estimator f(U,r],p) can be negative due to the stochastic esti- 
mation, the standard treatment is to absorb the sign into the observables, 
i.e. 

(Osgn(P))| P | 

{ ° )P ~ (sm(P)) lPl ■ (9) 

With the probability for the gauge link variable U and noise £ = (rj, p) 
written as P(U, £) oc P 1 (U)P 2 (U, £)P 3 (£) with 

Pi (17) oc e~ s ^ u) 
P2(U,0<x\f(U,0\ 

OC oo 

i=l fc=2 

the following two steps are needed to prove detailed balance |171 12U| . 

(a) Let Ti ([/, 17') be the ergodic Markov matrix satisfying detailed balance 
with respect to P u in other words Pi(U)Ti(U, U')dU = Pi{U')Ti{U' , U)dU' . 
Then the transition matrix 

T 12 (U,U') = T^U') min [l, ] (H) 

satisfies detailed balance with respect to the P\(JJ)P2{U,£) (with £ fixed). 

(b) The transition matrix 



T 23 &0 = P 3 (Omin[l, 



(12) 



satisfies detailed balance with respect to PziU, £)Ps(£) (with [/ fixed). 

From (a), (b) it follows that T12 and T23 keep the original distribution 
P([/, £) invariant and interleaving them will lead to an ergodic Markov process 
with the desired fixed point. 



4.1 Pade - Z 2 Estimator of Tr\nM with Unbiased Subtraction 

In Eq. (J7J, one needs to calculate Xi — Jjjln M (U)r)i in the stochastic series 
expansion of the fcrmion determinant. An efficient method is developed to 
calculate it ^H]. First of all, the logarithm is approximated using a Pade 
approximation, which after the partial fraction expansion, has the form 

N P 

\nM(U, k) r> R m (U) = b a I + J2b t (M(U, k) + c t iy 1 (13) 

i=i 

where Np is the order of the Pade approximation, and the constants bi and 
Cj are the Pade coefficients. In our implementation we have used an 11-th 
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order approximation whose coefficients are tabulated in 18;. The traces of 
InM are then estimated by evaluating bilinears of the form 7]'Rm(U)t], If the 
components of rj are chosen from the complex Z 2 group, then the contributions 
to the variance of these bilinears come only from off diagonal elements of 
Rm(U) [22 E3- I n this sense, Z 2 noise is optimal and has been applied to 
the calculation of nucleon matrix elements involving quark loops |23| . An 
effective method reducing the variance is to subtract off a linear combination 
of traceless operators from Rm(U) and to consider 

E[Tr R M {U),rj\ = rj* (R M (U) - a^) V . (14) 

Here the Oi are operators with Tr Oi — 0. Clearly since the Oi are traceless 
they do not bias the estimators. The on are constants that can be tuned to 
minimize the fluctuations in E[Tr Rm{U),tj]. 

With other types of noise such as Gaussian noise, the variance receives 
contributions from diagonal terms which one cannot subtract off. In this case, 
the unbiased subtraction scheme described here is ineffective. In practice, the 
Oi are constructed by taking traceless terms from the hopping parameter 
expansion for Af _1 (t/). It is shown for Wilson fermions on a 8 3 x 12 lattice 
at /3 = 5.6, these subtractions can reduce the noise coming from the terms 
(M(U) + Cj) _1 in equation ljT3|l by a factor as large as 37 for k = 0.150 with 
50 Z 2 noises Q2]. 

4.2 Implementation of the Noisy Monte Carlo Algorithm 

The noisy Monte Carlo algorithm has been implemented for the Wilson dy- 
namical fermion with pure gauge update (Kentucky Noisy Monte Carlo Al- 
gorithm) for an 8 4 lattice with /3 = 5.5 and k — 0.155 [201 ■ Several tricks are 
employed to reduce the fluctuations of the Tr In M estimate and increase the 
acceptance. These include shifting the TrlnM with a constant, A/3 shift |24| . 
and splitting the TrlnM with N 'fractional flavors'. After all these efforts, 
the results are shown to agree with those from the HMC simulation. However, 
the autocorrelation is very long and the acceptance rate is low. This has to 
do with the fact that Tr In M is an extensive quantity which is proportional 
to volume and the stochastic series expansion of e x converges for x < 6 for 
a sample with the size of ~ 10 3 — 10 4 . This is a stringent requirement which 
requires the fractional flavor number N > 15 for this lattice. This can be seen 
from the distribution of x — ^}2f{TrR* M (U) — X^Plag — Xq) /N in Fig. 3 which 
shows that taking N to be 15, 20, and 25, the largest x value is less than 6. 

As the volume increases, this fractional flavor needs to be larger to keep 
x smaller than 6 for a sample of the size ~ 10 3 — 10 4 . At the present volume 
(8 4 ), the autocorrelation is already much longer than that of HMC, it is going 
to be even less efficient for larger volumes. This is generic for the noisy Monte 
Carlo algorithm which scale with volume as V 2 , while HMC scales as V 5 ^ 4 . 
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x=E, ( Tr R M (U) - X, Plaq - xO, ) / N 
Fig. 3. Distributions of x for the three noisy simulations 



5 Hybrid Noisy Monte Carlo Algorithm — a New 
Proposal 

It is clear that the inefficiency of the noisy Monte Carlo algorithm for the 
fermion determinant is due to the large fluctuation of the TrlnM estima- 
tor from one gauge configuration to the next. We shall propose a combined 
Hybrid Monte Carlo (HMC) and Noisy Monte Carlo (NMC) to remove such 
fluctuations in the context of the finite density. 

With the baryon number projection discussed in Sec. |3| we can write the 
partition function for the finite baryon sector with B baryons as 



f + q an j. *trift^-uf fr 7 d0e~ l3B6 detM^M [9] 

J detM*M[9 = 0] 

(15) 

In this case, one can update the momentum p, the gauge link variable U 
and the pseudofermion field 4> via HMC and then interleave with NMC for 
updating the determinant ratio 

±- C d9e- aBe detM^ M[6] , 

detM^M[9 = 0] ' [ ' 

As described in Sec. 01 NMC involves two Metropolis accept/reject steps 
to update the ratio with the Pade - Z2 estimator of the Tr In difference of 
the determinants, i.e. Tr(1nM^M[9] - In M^M[0 = 0]). It is pointed out [2SJ 
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that for zero temperature, one can approximate the continuous integral over 
with a discrete sum incorporating triality zero projection ^2 ^] so that the 
partition function is a mixture of different Zg for different baryon number B. 
In other words, the approximation 



±- f" Me- i3BB detMiM[0] — ► -!- e^W? detM* M&^] (17) 

2tt J 3B n 3B N 

leads to the mixing of the baryon sector B with those of B ± Bm, B ± 2_Bjy.... 
If £> is small and TVb > B, then the partition will be dominated by Zb 
with small mixture from Zb±b n , Zb±2B n , ■■■■ For example, if we take B = 1 
and -Bat = 5, the discrete approximation gives an admixture of partition 
function with baryon number B = 1, 5, 11, —4, —9, .... At zero temperature, the 
partition function Zb behaves like e~ BmN tat , one expects that the mixing 
due to baryons other than B = 1 will be exponentially suppressed when 
m n N t a t > 1. 

Two points need to be stressed. First of all, it is crucial to project out the 
particle number in Eq. I|17H before the Metropolis accept /reject step in order 
to overcome the overlap problem. Secondly, given that the ratio R in Eq. 1|16|) 
is replaced with a discrete sum 

3-Bjv — 1 

X = _i_ <y e -i^Tr(\nM^M[2$Z)-\nM^M[0])^ ^ 



N 



k=0 



which involves the difference between the Tr In M[^f] and Tr In M^M[0], 
it takes out the fluctuation due to the gauge configuration which plagued the 
Kentucky Noisy Monte Carlo simulation in Sec. 21 Furthermore, the Tr In 
difference is expected to be 0(1) as seen from Fig. 1. If indeed is the case, it 
should lead to a better convergence of the stochastic series expansion in Eq. 
J7J and the algorithm scales with volume the same as HMC. 



5.1 Another Application 

Here we consider another possible application of the Hybrid Noisy Monte 
Carlo algorithm. HMC usually deals with two degenerate flavors. However, 
nature comes with 3 distinct light flavors - u,d and s. To consider u and d as 
separate flavors, one can perform HMC with two degenerate flavors at the d 
quark mass and then employ NMC to update the determinant ratio 

Ru _ DetM d DetM u =c TrqnM, l -lnM A ) Qgs 

Ud DetM\DetM d 

Since both the u and d masses are much smaller than Aqc d , TV(ln M u — In Md) 
should be small. If the Tr In difference is small enough (e.g. 0(1)) so that the 
acceptance rate is high, it could be a feasible algorithm for treating u and d 
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as distinct flavors so that realistic comparison with experiments can be done 
someday. It is shown recently that the Rational Hybrid Monte Carlo Algo- 
rithm (RHMC) |261 177] works efficiently for two flavor staggered fermions. It 
can be applied to single flavors for Wilson, domain wall, or overlap fermions 
at the cost of one pesudofermion for each flavor. We should point out that, 
in comparison, the Hybrid Noisy approach discussed here saves one pseud- 
ofcrmion, but at the cost of having to update the determinant ratio R u d in 
Eq. P]l. 

While we think that the Hybrid Noisy Monte Carlo algorithm proposed 
here might overcome the overlap and the low acceptance problems and the de- 
terminant detM[9] is real in this approach, the fact that the Fourier transform 
in Eq. Ijl7(l involves the baryon number B may still lead to a sign problem 
in the thermodynamic limit when B and V are large. However, as an initial 
attempt, we are more interested in finding out if the algorithm works for a 
small B such as 1 or 2 in a relatively small box. 
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